Chapter 4 Diversity analysis
genome_counts_filt <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
select(where(~ sum(.) > 0)) %>%
select(-AC85) %>%
rownames_to_column(., "genome")
genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt$genome)
table <- genome_counts_filt %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character,as.numeric) %>%
rownames_to_column(., "sample")
master_table <- sample_metadata %>%
mutate(sample=Tube_code) %>%
mutate(Tube_code=str_remove_all(Tube_code, "_a")) %>%
filter(Tube_code %in% table$sample) %>%
mutate(treatment = sub("^\\d+_", "", time_point))%>%
left_join(., table, by=join_by("Tube_code" =="sample"))
sample_metadata <- master_table %>%
select(1:13)
genome_counts_filt <- master_table %>%
select(12,14:323) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character,as.numeric) %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
dplyr::select(where(~ !all(. == 0)))%>%
rownames_to_column(., "genome")
genome_counts_filtering <- master_table %>%
select(12,14:323) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character,as.numeric) %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
dplyr::select(where(~ !all(. == 0)))
genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt$genome)
match_data(genome_counts_filtering,genome_tree)
genome_metadata <- genome_metadata %>%
filter(genome %in% genome_counts_filt$genome)4.1 Calculate diversities
4.1.1 Alpha diversity
# Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 0) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(richness = 1) %>%
rownames_to_column(var = "sample")
neutral <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(neutral = 1) %>%
rownames_to_column(var = "sample")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, tree = genome_tree) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(phylogenetic = 1) %>%
rownames_to_column(var = "sample")
genome_gifts1 <- genome_gifts[genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]
genome_counts_filt1 <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts1),]
genome_counts_filt1 <- genome_counts_filt1 %>%
remove_rownames() %>%
column_to_rownames(var = "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
rownames_to_column(., "genome")
dist <- genome_gifts1 %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
functional <- genome_counts_filt1 %>%
filter(genome %in% rownames(dist)) %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, dist = dist) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(functional = 1) %>%
rownames_to_column(var = "sample") %>%
mutate(functional = if_else(is.nan(functional), 1, functional))
# Merge all metrics
alpha_div <- richness %>%
full_join(neutral, by = join_by(sample == sample)) %>%
full_join(phylogenetic, by = join_by(sample == sample))%>%
full_join(functional, by = join_by(sample == sample))4.1.2 Beta diversity
beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 0)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1)
beta_q1p <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, tree = genome_tree)
genome_gifts1 <- genome_gifts[genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]
genome_counts_filt1 <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts1),]
dist <- genome_gifts1 %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
beta_q1f <- genome_counts_filt1 %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, dist = dist)4.2 Does the antibiotic administration change the microbial community?
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics") ) 4.2.1 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics") ) %>%
ggplot(aes(y = value, x = time_point, color=time_point, fill=time_point)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(values=treatment_colors)+
scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")Richness
Modelq0GLMMNB <- glmer.nb(neutral ~ time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(4.0048) ( log )
Formula: neutral ~ time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
223.7 229.6 -107.9 215.7 28
Scaled residuals:
Min 1Q Median 3Q Max
-1.4630 -0.6452 -0.1754 0.6566 2.6564
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0.07822 0.2797
Number of obs: 32, groups: individual, 18
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.9152 0.1558 18.72 < 2e-16 ***
time_point2_Antibiotics -0.9432 0.2144 -4.40 1.08e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
tm_pnt2_Ant -0.538
Neutral
Model_neutral <- lme(fixed = neutral ~ time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
225.017 230.6218 -108.5085
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 3.238645 7.593502
Fixed effects: neutral ~ time_point
Value Std.Error DF t-value p-value
(Intercept) 19.29561 2.000903 17 9.643449 0e+00
time_point2_Antibiotics -11.58649 2.716098 13 -4.265859 9e-04
Correlation:
(Intr)
time_point2_Antibiotics -0.631
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.86139926 -0.51381109 -0.09492589 0.56226694 1.96603543
Number of Observations: 32
Number of Groups: 18
Phylogenetic
Model_phylo <- lme(fixed = phylogenetic ~ time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
138.9014 144.5062 -65.45072
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 1.046827 1.691404
Fixed effects: phylogenetic ~ time_point
Value Std.Error DF t-value p-value
(Intercept) 5.461590 0.4814205 17 11.344740 0.0000
time_point2_Antibiotics -0.811302 0.6097318 13 -1.330588 0.2062
Correlation:
(Intr)
time_point2_Antibiotics -0.584
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.65239178 -0.50119332 -0.03450314 0.42332882 1.59503693
Number of Observations: 32
Number of Groups: 18
Functional
Model_funct <- lme(fixed = functional ~ time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_funct)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
-14.95179 -9.346998 11.47589
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 1.203543e-06 0.1504954
Fixed effects: functional ~ time_point
Value Std.Error DF t-value p-value
(Intercept) 1.4855915 0.03650050 17 40.70058 0.0000
time_point2_Antibiotics -0.0344778 0.05331239 13 -0.64671 0.5291
Correlation:
(Intr)
time_point2_Antibiotics -0.685
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.2430042 -0.5618554 0.2917987 0.6160260 2.0402200
Number of Observations: 32
Number of Groups: 18
4.2.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics")) %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.013523 0.0135230 2.0312 999 0.17
Residuals 30 0.199726 0.0066575
adonis2(richness ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.575038 | 0.1307754 | 4.513519 | 0.001 |
| Residual | 30 | 10.468804 | 0.8692246 | NA | NA |
| Total | 31 | 12.043842 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.03475 0.034752 2.5153 999 0.105
Residuals 30 0.41450 0.013817
adonis2(neutral ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.778963 | 0.1605587 | 5.738054 | 0.001 |
| Residual | 30 | 9.300869 | 0.8394413 | NA | NA |
| Total | 31 | 11.079832 | 1.0000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.31299 0.312989 16.486 999 0.002 **
Residuals 30 0.56954 0.018985
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(phylo ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.569966 | 0.2885492 | 12.16735 | 0.001 |
| Residual | 30 | 3.870929 | 0.7114508 | NA | NA |
| Total | 31 | 5.440895 | 1.0000000 | NA | NA |
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.18537 0.185373 4.9261 999 0.026 *
Residuals 30 1.12892 0.037631
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(func ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(func))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(func))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.121064 | 0.361186 | 16.96203 | 0.001 |
| Residual | 30 | 1.982777 | 0.638814 | NA | NA |
| Total | 31 | 3.103841 | 1.000000 | NA | NA |
4.2.3 Structural zeros
acclimation_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point == "1_Acclimation") %>%
dplyr::select(sample) %>%
pull()
antibiotics_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point == "2_Antibiotics") %>%
dplyr::select(sample) %>% pull()
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",subset_meta$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>%
mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>%
mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>%
mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>%
filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE) %>%
mutate(present = case_when(
all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
!all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
!all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Acclimation
structural_zeros %>%
filter(present=="acclimation") %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) # A tibble: 9 × 2
# Rowwise:
phylum sample_count
<chr> <int>
1 p__Bacillota_A 43
2 p__Bacteroidota 15
3 p__Bacillota 8
4 p__Pseudomonadota 8
5 p__Cyanobacteriota 3
6 p__Verrucomicrobiota 2
7 p__Bacillota_B 1
8 p__Bacillota_C 1
9 p__Fusobacteriota 1
Antibiotics
# A tibble: 4 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_7:bin_000003 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Proteus s__Proteus cibarius
2 LI1_2nd_7:bin_000001 p__Bacillota_A c__Clostridia o__Clostridiales f__Clostridiaceae g__Sarcina s__
3 AH1_2nd_7:bin_000055 p__Bacillota c__Bacilli o__Mycoplasmatales f__Mycoplasmoidaceae g__Ureaplasma s__
4 AH1_2nd_13:bin_000025 p__Bacillota_A c__Clostridia o__Christensenellales f__UBA3700 g__ s__
Community elements differences in structural zeros
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>
4.2.4 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics") )%>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "time_point",
rand_formula = "(1|individual)",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output$res %>%
dplyr::select(taxon, lfc_time_point2_Antibiotics, p_time_point2_Antibiotics) %>%
filter(p_time_point2_Antibiotics < 0.05) %>%
dplyr::arrange(p_time_point2_Antibiotics) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_time_point2_Antibiotics)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_point2_Antibiotics, y=forcats::fct_reorder(genome,lfc_time_point2_Antibiotics), fill=phylum)) + #forcats::fct_rev()
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
# coord_flip()+
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))Phyla of the significant MAGs
ancombc_rand_table_mag%>%
filter(lfc_time_point2_Antibiotics<0) %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) phylum sample_count
1 Bacteroidota 14
2 Bacillota_A 4
3 Bacillota 2
4 Campylobacterota 1
phylum genus species
1 Campylobacterota Helicobacter_J
2 Bacillota_A
3 Bacteroidota Bacteroides
4 Bacteroidota Odoribacter
5 Bacteroidota Bacteroides
6 Bacillota Coprobacillus
7 Bacillota
8 Bacteroidota Phocaeicola
9 Bacteroidota Bacteroides
10 Bacteroidota Phocaeicola
11 Bacteroidota Odoribacter
12 Bacteroidota Bacteroides
13 Bacteroidota Bacteroides
14 Bacteroidota Parabacteroides
15 Bacteroidota
16 Bacteroidota Parabacteroides
17 Bacillota_A
18 Bacteroidota Alistipes
19 Bacillota_A
20 Bacillota_A
21 Bacteroidota Odoribacter
ancombc_rand_table_mag%>%
filter(lfc_time_point2_Antibiotics<0) %>%
count(genus, name = "sample_count") %>%
arrange(desc(sample_count)) genus sample_count
1 6
2 Bacteroides 5
3 Odoribacter 3
4 Parabacteroides 2
5 Phocaeicola 2
6 Alistipes 1
7 Coprobacillus 1
8 Helicobacter_J 1
4.2.5 Differences in the functional capacity
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Antibiotics"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.2.5.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.331 0.0319
2 Antibiotics 0.244 0.134
MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94311, p-value = 0.09179
Wilcoxon rank sum exact test
data: value by treatment
W = 190, p-value = 0.01768
alternative hypothesis: true location shift is not equal to 0
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-Antibiotics)%>%
mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=treatment_colors) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")4.2.5.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.329 0.0319
2 Antibiotics 0.255 0.126
MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.95731, p-value = 0.2316
Wilcoxon rank sum exact test
data: value by treatment
W = 188, p-value = 0.02195
alternative hypothesis: true location shift is not equal to 0
function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))| Code_function | Acclimation | Antibiotics | Function | Difference | treatment |
|---|---|---|---|---|---|
| B06 | 0.68059020 | 0.47329940 | Organic anion biosynthesis | 0.20729080 | Acclimation |
| B02 | 0.59930820 | 0.41576970 | Amino acid biosynthesis | 0.18353850 | Acclimation |
| D02 | 0.38530780 | 0.20616790 | Polysaccharide degradation | 0.17913990 | Acclimation |
| S03 | 0.27145162 | 0.09403129 | Spore | 0.17742033 | Acclimation |
| B01 | 0.84159250 | 0.69078390 | Nucleic acid biosynthesis | 0.15080860 | Acclimation |
| B07 | 0.44505530 | 0.30423320 | Vitamin biosynthesis | 0.14082210 | Acclimation |
| B08 | 0.44618700 | 0.32094170 | Aromatic compound biosynthesis | 0.12524530 | Acclimation |
| D09 | 0.25519710 | 0.13446900 | Antibiotic degradation | 0.12072810 | Acclimation |
| B04 | 0.54481600 | 0.42941310 | SCFA biosynthesis | 0.11540290 | Acclimation |
| D03 | 0.29173710 | 0.20687250 | Sugar degradation | 0.08486460 | Acclimation |
| D05 | 0.28853330 | 0.22270070 | Amino acid degradation | 0.06583260 | Acclimation |
| B10 | 0.05960097 | 0.03635236 | Antibiotic biosynthesis | 0.02324861 | Acclimation |
4.3 Does antibiotic administration remove the differences found in the two populations?
4.3.1 Alpha diversity
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(time_point == "2_Antibiotics" )alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(time_point == "2_Antibiotics" ) %>%
ggplot(aes(y = value, x = Population, color=Population, fill=Population)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(values=treatment_colors)+
scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")4.3.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(time_point == "2_Antibiotics") %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(time_point == "2_Antibiotics")Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$Population) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.015319 0.0153186 6.8764 999 0.021 *
Residuals 21 0.046782 0.0022277
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(richness ~ Population,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.356644 | 0.1527052 | 3.784762 | 0.001 |
| Residual | 21 | 7.527428 | 0.8472948 | NA | NA |
| Total | 22 | 8.884073 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$Population) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.030536 0.0305358 3.8593 999 0.061 .
Residuals 21 0.166158 0.0079123
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(neutral ~ Population,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.785669 | 0.2085055 | 5.532084 | 0.001 |
| Residual | 21 | 6.778468 | 0.7914945 | NA | NA |
| Total | 22 | 8.564137 | 1.0000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$Population) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.012041 0.012041 0.9898 999 0.332
Residuals 21 0.255459 0.012165
adonis2(phylo ~ Population,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.8963254 | 0.1888758 | 4.889993 | 0.001 |
| Residual | 21 | 3.8492558 | 0.8111242 | NA | NA |
| Total | 22 | 4.7455811 | 1.0000000 | NA | NA |
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$Population) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.01804 0.018038 0.439 999 0.495
Residuals 21 0.86285 0.041088
adonis2(func ~ Population,
data = subset_meta %>% arrange(match(Tube_code,labels(func))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.0184644 | 0.0098378 | 0.2086464 | 0.704 |
| Residual | 21 | 1.8584184 | 0.9901622 | NA | NA |
| Total | 22 | 1.8768828 | 1.0000000 | NA | NA |
4.4 Are the microbial communities similar in both donor samples?
samples_to_keep <- sample_metadata %>%
filter(type =="Hot_control" & time_point %in% c( "3_Transplant1", "4_Transplant2"))%>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type =="Hot_control" & time_point %in% c( "3_Transplant1", "4_Transplant2"))alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(type =="Hot_control" & time_point %in% c( "3_Transplant1", "4_Transplant2")) %>%
ggplot(aes(y = value, x = time_point, color=time_point, fill=time_point)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(values=treatment_colors)+
scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000151 0.0001514 0.0181 999 0.905
Residuals 13 0.108675 0.0083596
adonis2(richness ~ time_point,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.06781753 | 0.02960432 | 0.3965972 | 0.979 |
| Residual | 13 | 2.22298047 | 0.97039568 | NA | NA |
| Total | 14 | 2.29079799 | 1.00000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00067 0.0006698 0.0523 999 0.801
Residuals 13 0.16645 0.0128035
adonis2(neutral ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.09261373 | 0.03872219 | 0.523666 | 0.804 |
| Residual | 13 | 2.29913452 | 0.96127781 | NA | NA |
| Total | 14 | 2.39174825 | 1.00000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00371 0.0037095 0.3076 999 0.635
Residuals 13 0.15676 0.0120588
adonis2(phylo ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.009763176 | 0.01920889 | 0.2546063 | 0.853 |
| Residual | 13 | 0.498500276 | 0.98079111 | NA | NA |
| Total | 14 | 0.508263452 | 1.00000000 | NA | NA |
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.013298 0.0132975 1.5346 999 0.217
Residuals 13 0.112644 0.0086649
adonis2(func ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(func))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | -0.003102642 | -0.01647707 | -0.2107297 | 0.812 |
| Residual | 13 | 0.191403206 | 1.01647707 | NA | NA |
| Total | 14 | 0.188300564 | 1.00000000 | NA | NA |
4.5 Does the donor sample maintain the microbial community found in acclimation?
4.5.1 Alpha diversity
sample_metadata <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Transplant1", "Transplant2") ~ "Transplant",
TRUE ~ treatment
))
samples_to_keep <- sample_metadata %>%
filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant"))alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(values=treatment_colors)+
scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")4.5.2 Beta diversity
Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.008956 0.0089559 1.5113 999 0.272
Residuals 22 0.130368 0.0059258
adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.2644898 | 0.06349573 | 1.491617 | 0.103 |
| Residual | 22 | 3.9009835 | 0.93650427 | NA | NA |
| Total | 23 | 4.1654732 | 1.00000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000794 0.0007940 0.0884 999 0.798
Residuals 22 0.197689 0.0089859
adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3163831 | 0.07591112 | 1.807234 | 0.059 |
| Residual | 22 | 3.8514266 | 0.92408888 | NA | NA |
| Total | 23 | 4.1678097 | 1.00000000 | NA | NA |
neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
# scale_color_manual(values = treatment_colors,labels=c("Cold_wet" = "Cold wet", "Hot_dry" = "Hot dry")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.002011 0.0020114 0.2582 999 0.693
Residuals 22 0.171393 0.0077906
adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.0409768 | 0.05601243 | 1.305392 | 0.286 |
| Residual | 22 | 0.6905894 | 0.94398757 | NA | NA |
| Total | 23 | 0.7315662 | 1.00000000 | NA | NA |
phylo %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000522 0.0005217 0.052 999 0.828
Residuals 22 0.220603 0.0100274
adonis2(func ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(func))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.08891682 | 0.2161114 | 6.065212 | 0.046 |
| Residual | 22 | 0.32252296 | 0.7838886 | NA | NA |
| Total | 23 | 0.41143978 | 1.0000000 | NA | NA |
func %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(values = treatment_colors,labels=c("Cold_wet" = "Cold wet", "Hot_dry" = "Hot dry")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)4.6 Is the donor sample microbiota different to recipient microbiota?
4.6.1 Alpha diversity
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant"))alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(values=treatment_colors)+
scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")4.6.2 Beta diversity
Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.12878 0.128776 13.013 999 0.001 ***
Residuals 20 0.19792 0.009896
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.76810 | 0.2746704 | 7.573672 | 0.001 |
| Residual | 20 | 4.66907 | 0.7253296 | NA | NA |
| Total | 21 | 6.43717 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.071939 0.071939 5.5113 999 0.038 *
Residuals 20 0.261057 0.013053
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 2.109391 | 0.3208221 | 9.447366 | 0.001 |
| Residual | 20 | 4.465565 | 0.6791779 | NA | NA |
| Total | 21 | 6.574956 | 1.0000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04724 0.047241 2.6056 999 0.111
Residuals 20 0.36261 0.018131
adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3764016 | 0.2391132 | 6.285119 | 0.001 |
| Residual | 20 | 1.1977547 | 0.7608868 | NA | NA |
| Total | 21 | 1.5741563 | 1.0000000 | NA | NA |
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000158 0.0001578 0.0131 999 0.9
Residuals 20 0.240836 0.0120418
adonis2(func ~ treatment,
data = subset_meta %>% arrange(match(Tube_code,labels(func))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.07897425 | 0.1809165 | 4.417536 | 0.084 |
| Residual | 20 | 0.35754884 | 0.8190835 | NA | NA |
| Total | 21 | 0.43652308 | 1.0000000 | NA | NA |
4.7 Does FMT change the microbial community over time?
4.7.1 Alpha diversity
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Treatment" & treatment %in% c("Post-FMT1","Post-FMT2") ) alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(type=="Treatment" & treatment %in% c("Post-FMT1", "Post-FMT2" )) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(values=treatment_colors)+
scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")Richness
Model_rich <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Model_rich)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(4.3157) ( log )
Formula: richness ~ treatment + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
171.5 174.8 -81.7 163.5 13
Scaled residuals:
Min 1Q Median 3Q Max
-1.72426 -0.71683 -0.08993 0.38000 1.84050
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0 0
Number of obs: 17, groups: individual, 9
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9416 0.1772 22.247 <2e-16 ***
treatmentPost-FMT2 0.4080 0.2420 1.686 0.0918 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
trtmnP-FMT2 -0.732
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
Neutral
Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
129.347 132.1792 -60.6735
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 3.521142 11.47471
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 24.81671 4.241888 8 5.850393 0.0004
treatmentPost-FMT2 15.06479 5.589803 7 2.695048 0.0309
Correlation:
(Intr)
treatmentPost-FMT2 -0.701
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.75054530 -0.55659466 -0.07097174 0.51610612 1.27053535
Number of Observations: 17
Number of Groups: 9
Phylogenetic
Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
56.51921 59.35141 -24.25961
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.7091595 0.8466237
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 4.192894 0.3867555 8 10.841199 0.0000
treatmentPost-FMT2 0.939941 0.4163443 7 2.257606 0.0585
Correlation:
(Intr)
treatmentPost-FMT2 -0.582
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.1453262 -0.5912067 -0.1412005 0.4845880 2.3194143
Number of Observations: 17
Number of Groups: 9
Functional
Model_funct <- lme(fixed = functional ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_funct)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
-32.8173 -29.9851 20.40865
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.0619285 0.03018711
Fixed effects: functional ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 1.4883350 0.02345764 8 63.44777 0.0000
treatmentPost-FMT2 0.0323233 0.01501285 7 2.15304 0.0683
Correlation:
(Intr)
treatmentPost-FMT2 -0.352
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.08983687 -0.66983137 -0.03874894 0.52553040 1.62217086
Number of Observations: 17
Number of Groups: 9
4.7.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Post-FMT1", "Post-FMT2")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Post-FMT1", "Post-FMT2"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.017295 0.0172950 2.8289 999 0.111
Residuals 15 0.091706 0.0061137
adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3571397 | 0.07959561 | 1.297184 | 0.00390625 |
| Residual | 15 | 4.1297869 | 0.92040439 | NA | NA |
| Total | 16 | 4.4869265 | 1.00000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.009753 0.0097526 0.7091 999 0.44
Residuals 15 0.206312 0.0137541
adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3157079 | 0.08117526 | 1.325203 | 0.00390625 |
| Residual | 15 | 3.5735051 | 0.91882474 | NA | NA |
| Total | 16 | 3.8892130 | 1.00000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.010945 0.010944 0.6593 999 0.529
Residuals 15 0.248985 0.016599
adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.06393539 | 0.09448624 | 1.565182 | 0.0703125 |
| Residual | 15 | 0.61272811 | 0.90551376 | NA | NA |
| Total | 16 | 0.67666350 | 1.00000000 | NA | NA |
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00162 0.0016197 0.1158 999 0.694
Residuals 15 0.20977 0.0139849
adonis2(func ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(func))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(func))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.02837271 | 0.07925317 | 1.291123 | 0.3203125 |
| Residual | 15 | 0.32962828 | 0.92074683 | NA | NA |
| Total | 16 | 0.35800099 | 1.00000000 | NA | NA |
4.8 Do FMT change the microbial community compared to antibiotics and acclimation?
4.8.1 Alpha diversity
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Treatment" & treatment %in% c("Antibiotics", "Acclimation","Post-FMT1") ) alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(type=="Treatment" & treatment %in% c( "Antibiotics","Acclimation", "Post-FMT1")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
# scale_color_manual(values=treatment_colors)+
# scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")Richness: Problems to calculate
Model_rich <- glmer.nb(richness ~ treatment + (1|individual), data = alpha_div_meta)
#summary(Model_rich)
emmeans(Model_rich, pairwise ~ treatment)Neutral
Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
#summary(Model_neutral)
emmeans(Model_neutral, pairwise ~ treatment)$emmeans
treatment emmean SE df lower.CL upper.CL
Acclimation 17.4 3.46 8 9.410 25.4
Antibiotics 8.5 3.91 8 -0.516 17.5
Post-FMT1 24.8 3.67 8 16.296 33.2
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Antibiotics 8.90 4.82 13 1.847 0.1935
Acclimation - (Post-FMT1) -7.36 4.62 13 -1.591 0.2840
Antibiotics - (Post-FMT1) -16.25 4.92 13 -3.301 0.0148
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates
Phylogenetic
Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
#summary(Model_phylo)
emmeans(Model_phylo, pairwise ~ treatment)$emmeans
treatment emmean SE df lower.CL upper.CL
Acclimation 5.56 0.514 8 4.38 6.75
Antibiotics 4.33 0.581 8 2.99 5.67
Post-FMT1 4.17 0.544 8 2.92 5.43
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Antibiotics 1.229 0.727 13 1.690 0.2456
Acclimation - (Post-FMT1) 1.387 0.698 13 1.986 0.1549
Antibiotics - (Post-FMT1) 0.158 0.744 13 0.212 0.9755
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates
Functional
Model_funct <- lme(fixed = functional ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
#summary(Model_funct)
emmeans(Model_funct, pairwise ~ treatment)$emmeans
treatment emmean SE df lower.CL upper.CL
Acclimation 1.49 0.0437 8 1.39 1.60
Antibiotics 1.47 0.0496 8 1.36 1.58
Post-FMT1 1.49 0.0464 8 1.38 1.60
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Antibiotics 0.02516 0.0661 13 0.380 0.9238
Acclimation - (Post-FMT1) 0.00591 0.0638 13 0.093 0.9953
Antibiotics - (Post-FMT1) -0.01924 0.0679 13 -0.283 0.9569
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates
4.8.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Antibiotics", "Post-FMT1")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Antibiotics", "Post-FMT1"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.015456 0.0077280 1.0306 999 0.4
Residuals 21 0.157471 0.0074986
adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.837434 | 0.201293 | 2.646248 | 0.001 |
| Residual | 21 | 7.290723 | 0.798707 | NA | NA |
| Total | 23 | 9.128157 | 1.000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.8277738 2.288600 0.1405032 0.001 0.003 *
2 Acclimation vs Post-FMT1 1 0.9551308 2.926054 0.1632291 0.001 0.003 *
3 Antibiotics vs Post-FMT1 1 0.9744543 2.741152 0.1741392 0.003 0.009 *
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.017834 0.0089172 0.5922 999 0.59
Residuals 21 0.316187 0.0150565
adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 2.356335 | 0.2697712 | 3.879056 | 0.001 |
| Residual | 21 | 6.378232 | 0.7302288 | NA | NA |
| Total | 23 | 8.734566 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.8812272 2.745592 0.1639591 0.001 0.003 *
2 Acclimation vs Post-FMT1 1 1.3127773 4.629826 0.2358567 0.001 0.003 *
3 Antibiotics vs Post-FMT1 1 1.3423457 4.351968 0.2508054 0.001 0.003 *
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.15470 0.077349 2.6532 999 0.098 .
Residuals 21 0.61221 0.029153
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.326050 | 0.3645754 | 6.024383 | 0.001 |
| Residual | 21 | 2.311195 | 0.6354246 | NA | NA |
| Total | 23 | 3.637246 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.6891376 5.128491 0.2681074 0.002 0.006 *
2 Acclimation vs Post-FMT1 1 0.2531943 3.273842 0.1791546 0.036 0.108
3 Antibiotics vs Post-FMT1 1 1.0996467 9.041595 0.4102060 0.003 0.009 *
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.23450 0.117249 4.3966 999 0.025 *
Residuals 21 0.56003 0.026668
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(func ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(func))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(func))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.636364 | 0.5726008 | 14.0672 | 0.002 |
| Residual | 21 | 1.221410 | 0.4273992 | NA | NA |
| Total | 23 | 2.857774 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 1.11238376 14.789477 0.5137112 0.001 0.003 *
2 Acclimation vs Post-FMT1 1 0.05911164 2.624236 0.1488993 0.137 0.411
3 Antibiotics vs Post-FMT1 1 1.36464611 16.864488 0.5647004 0.001 0.003 *
4.9 Is the gut microbiota similar to the donor after FMT?
4.9.1 Beta diversity
sample_metadata <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
TRUE ~ treatment
))
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.11342 0.113420 11.957 999 0.001 ***
Residuals 28 0.26559 0.009485
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.182046 | 0.154139 | 5.102365 | 0.001 |
| Residual | 28 | 6.486654 | 0.845861 | NA | NA |
| Total | 29 | 7.668700 | 1.000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04487 0.044865 3.2331 999 0.079 .
Residuals 28 0.38855 0.013877
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.146508 | 0.1607364 | 5.362581 | 0.001 |
| Residual | 28 | 5.986340 | 0.8392636 | NA | NA |
| Total | 29 | 7.132848 | 1.0000000 | NA | NA |
neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
# scale_color_manual(values = treatment_colors) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00008 0.0000787 0.005 999 0.937
Residuals 28 0.43669 0.0155959
adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.1296962 | 0.1018031 | 3.173567 | 0.002 |
| Residual | 28 | 1.1442941 | 0.8981969 | NA | NA |
| Total | 29 | 1.2739903 | 1.0000000 | NA | NA |
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00120 0.001202 0.1035 999 0.744
Residuals 28 0.32533 0.011619
adonis2(func ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(func))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(func))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.01544812 | 0.02751097 | 0.7920986 | 0.46 |
| Residual | 28 | 0.54607759 | 0.97248903 | NA | NA |
| Total | 29 | 0.56152571 | 1.00000000 | NA | NA |
4.9.2 Structural zeros
FMT_samples <- sample_metadata %>%
filter(type == "Treatment" & treatment == "FMT") %>%
dplyr::select(sample) %>%
pull()
Transplant_samples <- sample_metadata %>%
filter(type == "Treatment" & treatment =="Transplant") %>%
dplyr::select(sample) %>% pull()
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",subset_meta$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_FMT = all(c_across(all_of(FMT_samples)) == 0)) %>%
mutate(all_zeros_Transplant = all(c_across(all_of(Transplant_samples)) == 0)) %>%
mutate(average_FMT = mean(c_across(all_of(FMT_samples)), na.rm = TRUE)) %>%
mutate(average_Transplant = mean(c_across(all_of(Transplant_samples)), na.rm = TRUE)) %>%
filter(all_zeros_FMT == TRUE || all_zeros_Transplant==TRUE) %>%
mutate(present = case_when(
all_zeros_FMT & !all_zeros_Transplant ~ "Transplant",
!all_zeros_FMT & all_zeros_Transplant ~ "FMT",
!all_zeros_FMT & !all_zeros_Transplant ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "FMT", average_FMT, average_Transplant)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Structural zeros in each group
fmt <- structural_zeros %>%
filter(present=="FMT") %>%
count(phylum, name = "FMT") %>%
arrange(desc(FMT))
fmt_f <- structural_zeros %>%
filter(present=="FMT") %>%
count(family, name = "FMT") %>%
arrange(desc(FMT)) structural_zeros %>%
filter(present=="Transplant") %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt, by="phylum" ) %>%
mutate(across(everything(), ~ replace_na(., 0))) %>%
as.data.frame() %>%
summarise(across(where(is.numeric), ~ sum(.x, na.rm = TRUE))) Transplant FMT
1 52 55
Phyla to which the structural zeros belong in each group
structural_zeros %>%
filter(present=="Transplant") %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt, by="phylum" ) %>%
mutate(across(everything(), ~ replace_na(., 0))) %>%
tt()| phylum | Transplant | FMT |
|---|---|---|
| p__Bacillota_A | 20 | 21 |
| p__Bacillota | 16 | 10 |
| p__Pseudomonadota | 6 | 11 |
| p__Bacteroidota | 3 | 6 |
| p__Desulfobacterota | 2 | 2 |
| p__Bacillota_B | 1 | 0 |
| p__Chlamydiota | 1 | 0 |
| p__Cyanobacteriota | 1 | 0 |
| p__Fusobacteriota | 1 | 0 |
| p__Verrucomicrobiota | 1 | 2 |
| p__Actinomycetota | 0 | 1 |
| p__Bacillota_C | 0 | 1 |
| p__Campylobacterota | 0 | 1 |
Families to which the structural zeros belong in each group
structural_zeros %>%
filter(present=="Transplant") %>%
count(family, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt_f, by="family" ) %>%
mutate(across(everything(), ~ replace_na(., 0))) %>%
tt()| family | Transplant | FMT |
|---|---|---|
| f__Lachnospiraceae | 7 | 9 |
| f__Erysipelotrichaceae | 6 | 5 |
| f__UBA660 | 6 | 0 |
| f__Enterobacteriaceae | 5 | 2 |
| f__CAG-508 | 3 | 0 |
| f__Ruminococcaceae | 3 | 3 |
| f__Anaerovoracaceae | 2 | 0 |
| f__Coprobacillaceae | 2 | 0 |
| f__Desulfovibrionaceae | 2 | 2 |
| f__Oscillospiraceae | 2 | 1 |
| f__Tannerellaceae | 2 | 1 |
| f__UBA1242 | 2 | 0 |
| f__ | 1 | 3 |
| f__Akkermansiaceae | 1 | 0 |
| f__CAG-239 | 1 | 2 |
| f__Chlamydiaceae | 1 | 0 |
| f__Enterococcaceae | 1 | 3 |
| f__Fusobacteriaceae | 1 | 0 |
| f__Gastranaerophilaceae | 1 | 0 |
| f__Marinifilaceae | 1 | 0 |
| f__Mycoplasmoidaceae | 1 | 1 |
| f__Peptococcaceae | 1 | 0 |
| f__Anaerotignaceae | 0 | 2 |
| f__Bacteroidaceae | 0 | 2 |
| f__LL51 | 0 | 2 |
| f__UBA3700 | 0 | 2 |
| f__Acutalibacteraceae | 0 | 1 |
| f__Arcobacteraceae | 0 | 1 |
| f__CAG-274 | 0 | 1 |
| f__CALVMC01 | 0 | 1 |
| f__Devosiaceae | 0 | 1 |
| f__Mycobacteriaceae | 0 | 1 |
| f__Pumilibacteraceae | 0 | 1 |
| f__RUG11792 | 0 | 1 |
| f__Rhizobiaceae | 0 | 1 |
| f__Rikenellaceae | 0 | 1 |
| f__Sphingobacteriaceae | 0 | 1 |
| f__Streptococcaceae | 0 | 1 |
| f__UBA1997 | 0 | 1 |
| f__UBA3830 | 0 | 1 |
| f__Weeksellaceae | 0 | 1 |
4.9.2.1 Functional capacities of the structural zeros
#Aggregate bundle-level GIFTs into the compound level
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% structural_zeros$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
#Aggregate element-level GIFTs into the function level
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
#Aggregate function-level GIFTs into overall Biosynthesis, Degradation and Structural GIFTs
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
TRUE ~ treatment
))%>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
filter(genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=FMT-Transplant)%>%
mutate(group_color = ifelse(Difference <0, "Transplant","FMT")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=treatment_colors) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Elevation")4.9.3 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
TRUE ~ treatment
))%>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant")) %>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "treatment",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output$res %>%
dplyr::select(taxon, lfc_treatmentTransplant, p_treatmentTransplant) %>%
filter(p_treatmentTransplant < 0.05) %>%
dplyr::arrange(p_treatmentTransplant) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_treatmentTransplant)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()Differential abundance MAGs in each group
fmt_count <- ancombc_rand_table_mag %>%
filter(lfc_treatmentTransplant<0) %>%
count(phylum, name = "FMT") %>%
arrange(desc(FMT))
ancombc_rand_table_mag %>%
filter(lfc_treatmentTransplant>0) %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt_count, by="phylum")%>%
mutate(across(where(is.numeric), ~ replace_na(., 0))) phylum Transplant FMT
1 Bacteroidota 13 5
2 Bacillota_A 4 18
3 Bacillota 3 1
4 Campylobacterota 1 1
5 Desulfobacterota 1 2
6 Spirochaetota 1 0
7 Pseudomonadota 0 5
8 Cyanobacteriota 0 2
9 Bacillota_B 0 1
10 Bacillota_C 0 1
ancombc_rand_table_mag %>%
filter(lfc_treatmentTransplant>0) %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt_count, by="phylum")%>%
mutate(across(where(is.numeric), ~ replace_na(., 0))) %>%
as.data.frame() %>%
summarise(across(where(is.numeric), ~ sum(.x, na.rm = TRUE))) Transplant FMT
1 23 36
ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_treatmentTransplant, y=forcats::fct_reorder(genome,lfc_treatmentTransplant), fill=phylum)) + #forcats::fct_rev()
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
# coord_flip()+
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))4.9.4 Differences in the functional capacities
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
TRUE ~ treatment
))%>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.9.4.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 FMT 0.354 0.0255
2 Transplant 0.351 0.0426
MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94866, p-value = 0.1557
Wilcoxon rank sum exact test
data: value by treatment
W = 128, p-value = 0.4826
alternative hypothesis: true location shift is not equal to 0
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=FMT-Transplant)%>%
mutate(group_color = ifelse(Difference <0, "Transplant","FMT")) means_gift <- element_gift_filt %>%
select(-treatment) %>%
pivot_longer(!sample, names_to = "Elements", values_to = "abundance") %>%
left_join(sample_metadata, by=join_by(sample==sample)) %>%
group_by(treatment, Elements) %>%
summarise(mean=mean(abundance))
log_fold <- means_gift %>%
group_by(Elements) %>%
summarise(
logfc_fmt_transplant = log2(mean[treatment == "FMT"] / mean[treatment == "Transplant"])
) %>%
left_join(., difference_table, by="Elements")difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=treatment_colors) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")log_fold %>%
filter(Elements!="D0611") %>%
ggplot(aes(x=forcats::fct_reorder(Function,logfc_fmt_transplant), y=logfc_fmt_transplant, fill=group_color)) +
geom_col() +
scale_fill_manual(values=treatment_colors) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Log_fold")+
labs(fill="Treatment")4.9.4.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 FMT 0.349 0.0221
2 Transplant 0.354 0.0374
MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.87186, p-value = 0.001843
Wilcoxon rank sum exact test
data: value by treatment
W = 121, p-value = 0.6801
alternative hypothesis: true location shift is not equal to 0
function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))